2.1_multiome_mouse_mammary_gland
Anna Schwager
2025-02-28
Load the global variables
here::i_am("scripts/2.1_multiome_mouse_mammary_gland.Rmd")
mainDir = here::here()
source(knitr::purl(file.path(mainDir,"scripts/global_variables.Rmd"), quiet=TRUE))
source(knitr::purl(file.path(mainDir, "scripts/functions.Rmd"), quiet=TRUE))
set.seed(123)1 Preparation and object creation
1.1 Setting directories and input files
inputDir_local = file.path(inputDir, "mm10", "mouse_mammary_gland", "one_cell_multiome")
outputDir_local = file.path(outputDir, "2.1_multiome_mouse_mammary_gland") ; if(!file.exists(outputDir_local)){dir.create(outputDir_local)}
outputDir_objects = file.path(outputDir_local, "objects") ; if(!file.exists(outputDir_objects)){dir.create(outputDir_objects)}
outputDir_plots = file.path(outputDir_local, "plots") ; if(!file.exists(outputDir_plots)){dir.create(outputDir_plots)}
outputDir_tables = file.path(outputDir_local, "tables") ; if(!file.exists(outputDir_tables)){dir.create(outputDir_tables)}
ld <- list.dirs(inputDir_local)
fragpaths <- list.files(ld[grepl("/.*fragmentFiles", ld)], full.names = TRUE)
fragpaths_k4 <- fragpaths[grepl(".*/h3k4me1/.*fragmentFiles/.*tsv.gz$", fragpaths)]
fragpaths_k27 <- fragpaths[grepl(".*/h3k27me3/.*fragmentFiles/.*tsv.gz$", fragpaths)]
fragpaths <- list(fragpaths_k27, fragpaths_k4)
names(fragpaths) <- c("h3k27me3", "h3k4me1")
rnapaths <- ld[grepl("/10XlikeMatrix_umi", ld)]
rm(ld)1.2 Loading annotation
consensus_peaks <- list("h3k27me3" = toGRanges(file.path(annotDir, "CREneg_peaks_h3k27me3.bed"), format="BED", header=FALSE),
"h3k4me1" = toGRanges(file.path(annotDir, "CREneg_peaks_h3k4me1.bed"), format="BED", header=FALSE))1.3 Loading the bigwigs and tables
ld <- list.dirs(inputDir)
bw_paths <- list.files(ld[grepl(".*/mouse_mammary_gland/bigwig.*", ld)], full.names = TRUE)
names(bw_paths) <- sub(".*D1888_CRE3-Mice8724_", "", bw_paths)
list_rna_bw <- list(rna_all = bw_paths[["all_rna_pseudobulk.bw"]],
rna_luminal = bw_paths[["cluster_1_rna_luminal.bw"]],
rna_basal = bw_paths[["cluster_0_rna_basal.bw"]],
rna_luminal_sc = bw_paths[["T0302_luminal_cell.bw"]],
rna_basal_sc = bw_paths[["T0290_basal_cell.bw"]])
rm(ld)1.4 Loading the FACS data
ld <- list.dirs(inputDir_local)
facs_path <- list.files(ld[grepl("/facs", ld)], full.names = TRUE)
facs_data <- read.csv(facs_path)
rm(ld)
row.names(facs_data) <- facs_data$CB
facs_data <- facs_data[, c(1,2)]
colnames(facs_data) <- c("CD24", "CD49f") 1.5 Creating the object
## loading RNA data
rna.data <- Read10X(data.dir = rnapaths)
rna_seurat <- CreateSeuratObject(counts = rna.data,
min.cells = 1,
min.features = -1,
project = "mouse_mammary_gland")
## adding FACS annotation
rna_seurat <- AddMetaData(
object = rna_seurat,
metadata = facs_data
)
## creating the multiome seurat objects
seurat_list <- list()
for (i in names(fragpaths)){
message(paste0("### Making fragment object for ", i))
total_counts <- CountFragments(fragpaths[[i]])
barcodes <- total_counts$CB
frags <- CreateFragmentObject(path = fragpaths[[i]], cells = barcodes)
message(paste0("### Making 50k bin matrix for ", i))
bin50k_kmatrix = GenomeBinMatrix(
frags,
genome = mm10,
cells = NULL,
binsize = 50000,
process_n = 10000,
sep = c(":", "_"),
verbose = TRUE)
message(paste0("### Making peak matrix for ", i))
peak_matrix = FeatureMatrix(
frags,
features = consensus_peaks[[i]],
cells = NULL,
sep = c("-", "-"),
verbose = TRUE)
message(paste0("### Creating chromatin assay for ", i))
chrom_assay <- CreateChromatinAssay(
counts = bin50k_kmatrix,
sep = c(":", "_"),
genome = "mm10",
fragments = fragpaths[[i]],
min.cells = 1,
min.features = -1)
message(paste0("### Creating Seurat object for ", i))
seurat <- CreateSeuratObject(
counts = chrom_assay,
assay = "bin_50k")
message(paste0("### Adding peak assay for ", i))
seurat[["peaks"]] <- CreateChromatinAssay(
counts = peak_matrix, genome = "mm10")
message(paste0("### Adding RNA assay for ", i))
seurat <- AddMetaData(seurat, rna_seurat@meta.data)
seurat[["RNA"]] <- rna_seurat@assays[["RNA"]]
Annotation(seurat) <- annotations_mm10
seurat <- AddMetaData(seurat, CountFragments(fragpaths[[i]]))
seurat <- FRiP(object = seurat, assay = "peaks", total.fragments = "frequency_count")
seurat@meta.data[["orig.ident"]] <- i
seurat_list[[i]] <- seurat
rm(seurat)
rm(frags)
rm(total_counts)
rm(barcodes)
rm(fragments_per_cell)
rm(bin50k_matrix)
rm(peak_matrix)
rm(chromatin_assay)
}
saveRDS(seurat_list, file.path(outputDir_objects, "seurat_list.rds"))
## merging h3k4me1 and h3k27me3 data into one object
seurat <- merge(seurat_list[[1]], seurat_list[[2]])
seurat <- JoinLayers(seurat, assay = "RNA")
saveRDS(seurat, file.path(outputDir_objects, "seurat_cre_step1.rds"))1.6 Reload the object - step1
seurat <- readRDS(file.path(outputDir_objects, "seurat_cre_step1.rds"))
seurat_list <- readRDS(file.path(outputDir_objects, "seurat_list.rds"))2 FACS annotation
2.1 Gating
facs_data <- facs_data %>%
mutate(phenotype = case_when(
(CD24 > 100000) | (CD24 > 25000 & CD49f < 8900) ~ "luminal",
TRUE ~ "basal"
))
ggplot(facs_data, aes(x = CD49f, y = CD24, color = phenotype)) +
geom_point(size = 3) +
labs(x = "CD49f", y = "CD24", title = "CD24 vs CD49f by Phenotype") +
theme_minimal() +
scale_color_manual(values = c("luminal" = mypal[5], "basal" = mypal[4]))
ggplot(facs_data, aes(x = log(CD49f), y = log(CD24), color = phenotype)) +
geom_point(size = 3) +
labs(x = "log(CD49f)", y = "log(CD24)", title = "CD24 vs CD49f by Phenotype log") +
theme_minimal() +
scale_color_manual(values = c("luminal" = mypal[5], "basal" = mypal[4]))
2.2 Adding FACS data as Seurat metadata
##the data will appear as phenotype column
seurat <- AddMetaData(
object = seurat,
metadata = facs_data
)3 QC
3.1 Weighted histograms
for (s in seurat_list){
p1 <- plot_weighted_hist(s) + ggtitle(print(s@meta.data[["orig.ident"]][1]), " DNA fragments")
p2 <- plot_weighted_hist(s, assay = "RNA") + ggtitle(print(s@meta.data[["orig.ident"]][1]), " RNA reads")
p3 <- plot_weighted_hist(s, assay = "RNA_features") + ggtitle(print(s@meta.data[["orig.ident"]][1]), " RNA features")
print(p1)
print(p2)
print(p3)
rm(s)
}## [1] "h3k27me3"
## [1] "h3k27me3"
## [1] "h3k27me3"## [1] "h3k4me1"
## [1] "h3k4me1"
## [1] "h3k4me1"
### Filtering
min_reads_chromatin = 400
min_reads_rna = 1000
seurat@meta.data[["filtering"]] <- ifelse(seurat@meta.data$frequency_count > min_reads_chromatin &
seurat@meta.data$nCount_RNA > min_reads_rna,
'pass', 'fail')3.2 FrIP and N fragments plots
p_frip <- VlnPlot(subset(seurat, subset = filtering == 'pass'), c("FRiP"),
group.by = "orig.ident", split.by = NULL, pt.size = 0) +
labs(title = "FRiP CRE mice") +
stat_summary(fun.y = median, geom='point', size = 2, colour = "black") +
theme(legend.position = "none") +
scale_fill_manual(values = c(mypal[3], mypal[2]))
print(p_frip)
p_count <- VlnPlot(subset(seurat, subset = filtering == 'pass'), c("frequency_count"),
group.by = "orig.ident", split.by = NULL, pt.size = 0, y.max = 12000) +
labs(title = "Unique fragments per cell CRE") +
stat_summary(fun.y = median, geom='point', size = 2, colour = "black") +
theme(legend.position = "none") +
scale_fill_manual(values = c(mypal[3], mypal[2]))
print(p_count)
p_count_rna <- VlnPlot(subset(seurat, subset = filtering == 'pass'), c("nCount_RNA"),
group.by = "orig.ident", split.by = NULL, pt.size = 0) +
labs(title = "Unique reads per cell RNA CRE") +
stat_summary(fun.y = median, geom='point', size = 2, colour = "black") +
theme(legend.position = "none") +
scale_fill_manual(values = c(mypal[3], mypal[2]))
print(p_count_rna)
p_genes_rna <- VlnPlot(subset(seurat, subset = filtering == 'pass'), c("nFeature_RNA"),
group.by = "orig.ident", split.by = NULL, pt.size = 0) +
labs(title = "Unique genes per cell CRE") +
stat_summary(fun.y = median, geom='point', size = 2, colour = "black") +
theme(legend.position = "none") +
scale_fill_manual(values = c(mypal[3], mypal[2]))
print(p_genes_rna)
ggsave(file.path(outputDir_plots, "frip_mouse_mammary_gland.pdf"),
plot = p_frip,
device = "pdf",
units = "mm",
width = 100,
height = 130)
ggsave(file.path(outputDir_plots, "frip_mouse_mammary_gland.pdf"),
plot = p_count,
device = "pdf",
units = "mm",
width = 100,
height = 130)
ggsave(file.path(outputDir_plots, "frip_mouse_mammary_gland.pdf"),
plot = p_genes_rna,
device = "pdf",
units = "mm",
width = 100,
height = 130)
ggsave(file.path(outputDir_plots, "frip_mouse_mammary_gland.pdf"),
plot = p_count_rna,
device = "pdf",
units = "mm",
width = 100,
height = 130)3.3 DNA to RNA reads
seurat <- subset(seurat, subset = filtering == "pass")
frags_to_rna_reads <- ggplot(seurat@meta.data, aes(x = frequency_count, y = nCount_RNA, color = orig.ident)) +
geom_point(alpha = 1, stroke = 0, size = 2) +
labs(x = "N unique DNA fragments",
y = "N unique RNA reads",
title = "N unique DNA fragments vs. N unique RNA reads") +
theme_classic() +
scale_color_manual(values = c(mypal[3], mypal[2]))
frags_to_genes <- ggplot(seurat@meta.data, aes(x = frequency_count, y = nFeature_RNA, color = orig.ident)) +
geom_point(alpha = 1, stroke = 0, size = 2) +
labs(x = "N unique DNA fragments",
y = "N unique genes",
title = "N unique genes vs. N unique RNA reads") +
theme_classic() +
scale_color_manual(values = c(mypal[3], mypal[2]))
genes_to_rna_reads <- ggplot(seurat@meta.data, aes(x = nFeature_RNA, y = nCount_RNA, color = orig.ident)) +
geom_point(alpha = 1, stroke = 0, size = 2) +
labs(x = "N unique genes",
y = "N unique RNA reads",
title = "N unique genes vs. N unique RNA reads") +
theme_classic() +
scale_color_manual(values = c(mypal[3], mypal[2]))
frags_to_rna_reads
frags_to_genes
genes_to_rna_reads
ggsave(file.path(outputDir_plots, "frags_to_rna_reads.pdf"),
plot = frags_to_rna_reads,
device = "pdf",
units = "mm",
width = 120,
height = 80)
ggsave(file.path(outputDir_plots, "frags_to_genes.pdf"),
plot = frags_to_genes,
device = "pdf",
units = "mm",
width = 120,
height = 80)
ggsave(file.path(outputDir_plots, "genes_to_rna_reads.pdf"),
plot = genes_to_rna_reads,
device = "pdf",
units = "mm",
width = 120,
height = 80)
4 Clustering
4.1 Clustering by RNA
seurat <- subset(seurat, subset = filtering == 'pass')
DefaultAssay(seurat) <- "RNA"
seurat <- SCTransform(seurat)
seurat <- RunPCA(seurat)
seurat <- RunUMAP(seurat, dims = 1:30, verbose = FALSE, seed.use = 42)
seurat <- FindNeighbors(seurat, dims = 1:30, verbose = FALSE)
seurat <- FindClusters(seurat, resolution = 0.2)
DimPlot(seurat, cols = c(mypal[11], mypal[10]))
saveRDS(seurat, file.path(outputDir_objects, "seurat_cre_processed_step2.rds"))seurat <- readRDS(file.path(outputDir_objects, "seurat_cre_processed_step2.rds"))4.2 Clustering by epigenomic mark
Separating h3k27me3 and h3k4me1 data
cre_list <- list("h3k4me1" = subset(seurat, subset = orig.ident == "h3k4me1"),
"h3k27me3" = subset(seurat, subset = orig.ident == "h3k27me3"))Running normalization, dim reduction and checking the principal components.
for (mark in names(cre_list)){
DefaultAssay(cre_list[[mark]]) <- 'bin_50k'
cre_list[[mark]] <- RunTFIDF(cre_list[[mark]], assay = 'bin_50k')
cre_list[[mark]] <- FindTopFeatures(cre_list[[mark]], min.cutoff = 'q0')
cre_list[[mark]] <- RunSVD(cre_list[[mark]])
p <- DepthCor(cre_list[[mark]])
print(p)
}
Component 1 should not be used.
for (mark in names(cre_list)){
cre_list[[mark]] <- RunUMAP(object = cre_list[[mark]],
reduction = 'lsi',
dims = 2:30)
cre_list[[mark]] <- FindNeighbors(object = cre_list[[mark]],
reduction = 'lsi',
dims = 2:30)
cre_list[[mark]] <- FindClusters(object = cre_list[[mark]],
algorithm = 3,
resolution = 0.8)
#p <- DimPlot(cre_list[[mark]])
#print(p)
}
saveRDS(cre_list, file.path(outputDir_objects, "seurat_cre_list_by_mark_processed_step2.rds"))cre_list <- readRDS(file.path(outputDir_objects, "seurat_cre_list_by_mark_processed_step2.rds"))Quick UMAPs to see sequencing-FACS correspondence
DimPlot(object = seurat, label = TRUE, cols = c(mypal[11], mypal[10])) + labs(title = "RNA clusters")
DimPlot(object = seurat, label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "RNA clusters, FACS annotation")
DimPlot(object = cre_list[["h3k27me3"]], label = TRUE, cols = c(mypal[11], mypal[10])) + labs(title = "H3K27me3 clusters")
DimPlot(object = cre_list[["h3k27me3"]], label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "H3K27me3 clusters, FACS annotation")
DimPlot(object = cre_list[["h3k4me1"]], label = TRUE, cols = c(mypal[11], mypal[10])) + labs(title = "H3K4me1 clusters")
DimPlot(object = cre_list[["h3k4me1"]], label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "H3K4me1 clusters, FACS annotation")
In RNA and K27: cluster 0 = basal-like, cluster 1 = luminal-like. In K4
it is vice versa.
Adding this annotation to meta_data:
tmp_vector <- cre_list[["h3k27me3"]]@meta.data[["seurat_clusters"]]
annot_k27 <- ifelse(tmp_vector == 1, "luminal_like_k27", "basal_like_k27")
cre_list[["h3k27me3"]]@meta.data[["annot_k27"]] <- annot_k27
tmp_vector2 <- cre_list[["h3k27me3"]]@meta.data[["SCT_snn_res.0.2"]]
annot_rna <- ifelse(tmp_vector2 == 1, "luminal_rna", "basal_rna")
cre_list[["h3k27me3"]]@meta.data[["annot_rna"]] <- annot_rna
tmp_vector <- cre_list[["h3k4me1"]]@meta.data[["seurat_clusters"]]
annot_k4 <- ifelse(tmp_vector == 1, "basal_like_k4", "luminal_like_k4")
cre_list[["h3k4me1"]]@meta.data[["annot_k4"]] <- annot_k4
tmp_vector2 <- cre_list[["h3k4me1"]]@meta.data[["SCT_snn_res.0.2"]]
annot_rna <- ifelse(tmp_vector2 == 1, "luminal_rna", "basal_rna")
cre_list[["h3k4me1"]]@meta.data[["annot_rna"]] <- annot_rna
cre_list[["h3k4me1"]]<- SetIdent(cre_list[["h3k4me1"]], value = "annot_k4")
cre_list[["h3k27me3"]]<- SetIdent(cre_list[["h3k27me3"]], value = "annot_k27")
tmp_vector3 <- seurat@meta.data[["SCT_snn_res.0.2"]]
annot_rna <- ifelse(tmp_vector3 == 1, "luminal_rna", "basal_rna")
seurat@meta.data[["annot_rna"]] <- annot_rna
seurat <- SetIdent(seurat, value = "annot_rna")
saveRDS(cre_list, file.path(outputDir_objects, "seurat_cre_list_by_mark_processed_step3.rds"))
saveRDS(seurat, file.path(outputDir_objects, "seurat_cre_processed_step3.rds"))Reload objects - step3
cre_list <- readRDS(file.path(outputDir_objects, "seurat_cre_list_by_mark_processed_step3.rds"))
seurat <- readRDS(file.path(outputDir_objects, "seurat_cre_processed_step3.rds"))4.3 UMAP plots
4.3.1 UMAPs colored by different annotations
cre_umap_rna <- DimPlot(object = seurat, label = TRUE, cols = c(mypal[11], mypal[10]), group.by = "annot_rna") + labs(title = "RNA clusters")
cre_umap_rna_facs <- DimPlot(object = seurat, label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "RNA clusters, FACS annotation")
cre_umap_k4 <- DimPlot(object = cre_list[["h3k4me1"]], label = TRUE, cols = c(mypal[11], mypal[10]), group.by = "annot_k4") + labs(title = "H3K4me1 clusters")
cre_umap_k4_facs <- DimPlot(object = cre_list[["h3k4me1"]], label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "H3K4me1 clusters, FACS annotation")
cre_umap_k27 <- DimPlot(object = cre_list[["h3k27me3"]], label = TRUE, cols = c(mypal[11], mypal[10]), group.by = "annot_k27") + labs(title = "H3K27me3 clusters")
cre_umap_k27_facs <- DimPlot(object = cre_list[["h3k27me3"]], label = TRUE, cols = c(mypal[4], mypal[5]), group.by = "phenotype") +
labs(title = "H3K27me3 clusters, FACS annotation")
cre_umap_rna
cre_umap_rna_facs
cre_umap_k4
cre_umap_k4_facs
cre_umap_k27
cre_umap_k27_facs
ggsave(file.path(outputDir_plots, "cre_umap_rna.pdf"),
plot = cre_umap_rna,
device = "pdf",
units = "mm",
width = 120,
height = 100)
ggsave(file.path(outputDir_plots, "cre_umap_rna_facs.pdf"),
plot = cre_umap_rna_facs,
device = "pdf",
units = "mm",
width = 120,
height = 100)
ggsave(file.path(outputDir_plots, "cre_umap_k4.pdf"),
plot = cre_umap_k4,
device = "pdf",
units = "mm",
width = 100,
height = 100)
ggsave(file.path(outputDir_plots, "cre_umap_k4_facs.pdf"),
plot = cre_umap_k4_facs,
device = "pdf",
units = "mm",
width = 80,
height = 100)
ggsave(file.path(outputDir_plots, "cre_umap_k27.pdf"),
plot = cre_umap_k27,
device = "pdf",
units = "mm",
width = 100,
height = 100)
ggsave(file.path(outputDir_plots, "cre_umap_k27_facs.pdf"),
plot = cre_umap_k27_facs,
device = "pdf",
units = "mm",
width = 80,
height = 100)4.3.2 Connected UMAP
# Extract UMAP coordinates from each representation
df1 <- cre_list[["h3k4me1"]]@reductions[["umap"]]@cell.embeddings %>% as.data.frame()
df2 <- seurat@reductions[["umap"]]@cell.embeddings %>% as.data.frame()
df3 <- cre_list[["h3k27me3"]]@reductions[["umap"]]@cell.embeddings %>% as.data.frame()
# Extract information about cell type or cluster
df1$cluster <- as.vector(cre_list[["h3k4me1"]]@meta.data[["annot_k4"]])
df2$cluster <- as.vector(seurat@meta.data[["annot_rna"]])
df3$cluster <- as.vector(cre_list[["h3k27me3"]]@meta.data[["annot_k27"]])
# Add the info about the mark
df1$mark <- "H3K4me1"
df2$mark <- "RNA"
df3$mark <- "H3K27me3"
# Shift umap_1 for each dataset
df1$umap_1 <- df1$umap_1 - 15 # Shift H3K4me1 left
df3$umap_1 <- df3$umap_1 + 20 # Shift H3K27me3 right
# Create a list and bind dataframes into one
dat.umap.lst <- list(df1, df2, df3)
dat.umap.long <- dat.umap.lst %>% bind_rows()
dat.umap.long$cell <- c(rownames(df1), rownames(df2), rownames(df3))
# Define color palette for clusters
cbPalette.ctype <- rep(c(mypal[10], mypal[11]), 3)
names(cbPalette.ctype) <- unique(dat.umap.long$cluster)
connected_umap_cre <- ggplot(dat.umap.long, aes(x = umap_1, y = umap_2, color = cluster, group = cell)) +
geom_point(alpha = 1) +
geom_path(alpha = 0.1) +
scale_color_manual(values = cbPalette.ctype) +
theme_classic() + ggtitle("H3K4me1 -> RNA -> H3K27me3")
connected_umap_cre 
ggsave(file.path(outputDir_plots, "connected_umap_cre_larger.pdf"),
plot = connected_umap_cre,
device = "pdf",
units = "mm",
width = 260,
height = 100)
ggsave(file.path(outputDir_plots, "connected_umap_cre_smaller.pdf"),
plot = connected_umap_cre,
device = "pdf",
units = "mm",
width = 220,
height = 80)5 Differential analysis
5.1 RNA cluster markers
rna_markers <- FindMarkers(seurat, assay = "SCT", ident.1 = "luminal_rna")
DotPlot(seurat,
features = rownames(rna_markers)[1:20],
cols = c("#EDE9E7", "#7e6148")) +
coord_flip() + labs(title = "Luminal vs basal markers, RNA") 
write.csv(rna_markers, file.path(outputDir_tables, "cre_rna_cluster_markers.csv"))5.2 Differentially “acessible” bins and peaks - basal vs luminal
DefaultAssay(cre_list[["h3k4me1"]]) <- 'bin_50k'
DefaultAssay(cre_list[["h3k27me3"]]) <- 'bin_50k'
da_bins_k4 <- FindMarkers(
object = cre_list[["h3k4me1"]],
ident.1 = "luminal_like_k4",
ident.2 = "basal_like_k4",
test.use = 'wilcox',
min.pct = 0.1
)
da_bins_k4$bin <- rownames(da_bins_k4)
da_bins_k4$query_region <- rownames(da_bins_k4)
da_bins_k27 <- FindMarkers(
object = cre_list[["h3k27me3"]],
ident.1 = "luminal_like_k27",
ident.2 = "basal_like_k27",
test.use = 'wilcox',
min.pct = 0.1
)
da_bins_k27$bin <- rownames(da_bins_k27)
da_bins_k27$query_region <- rownames(da_bins_k27)
DefaultAssay(cre_list[["h3k4me1"]]) <- 'peaks'
DefaultAssay(cre_list[["h3k27me3"]]) <- 'peaks'
da_peaks_k4 <- FindMarkers(
object = cre_list[["h3k4me1"]],
ident.1 = "luminal_like_k4",
ident.2 = "basal_like_k4",
test.use = 'wilcox',
min.pct = 0.1
)
da_peaks_k4$query_region <- rownames(da_peaks_k4) #to do inner_join with the closest feature
colnames(da_peaks_k4) <- c("p_val", "avg_log2FC", "pct.luminal_like", "pct.basal_like", "p_val_adj", "query_region")
da_peaks_k27 <- FindMarkers(
object = cre_list[["h3k27me3"]],
ident.1 = "luminal_like_k27",
ident.2 = "basal_like_k27",
test.use = 'wilcox',
min.pct = 0.1
)
da_peaks_k27$query_region <- rownames(da_peaks_k27)
colnames(da_peaks_k27) <- c("p_val", "avg_log2FC", "pct.luminal_like", "pct.basal_like", "p_val_adj", "query_region")Finding closest features.
features_da_k27 <- ClosestFeature(cre_list[["h3k27me3"]], regions = rownames(da_peaks_k27), annotation = annotations_mm10)
features_da_k4 <- ClosestFeature(cre_list[["h3k4me1"]], regions = rownames(da_peaks_k4), annotation = annotations_mm10)
da_peaks_k27 <- inner_join(da_peaks_k27, features_da_k27, by = "query_region")
da_peaks_k4 <- inner_join(da_peaks_k4, features_da_k4, by = "query_region")
features_da_k27_bins <- ClosestFeature(cre_list[["h3k27me3"]], regions = rownames(da_bins_k27), annotation = annotations_mm10)
features_da_k4_bins <- ClosestFeature(cre_list[["h3k4me1"]], regions = rownames(da_bins_k4), annotation = annotations_mm10)
da_bins_k27 <- inner_join(da_bins_k27, features_da_k27_bins, by = "query_region")
da_bins_k4 <- inner_join(da_bins_k4, features_da_k4_bins, by = "query_region")
write.csv(da_peaks_k27, file.path(outputDir_tables, "cre_diff_peaks_with_closest_genes_k27.csv"))
write.csv(da_peaks_k4, file.path(outputDir_tables, "cre_diff_peaks_with_closest_genes_k4.csv"))
write.csv(da_bins_k27, file.path(outputDir_tables, "cre_diff_bins_with_closest_genes_k27.csv"))
write.csv(da_bins_k4, file.path(outputDir_tables, "cre_diff_bins_with_closest_genes_k4.csv"))intersect(row.names(rna_markers[rna_markers$p_val_adj < 0.05, ]),
da_bins_k27[da_bins_k27$p_val < 0.05, ]$gene_name)
#Top 3: "Cxcl14" "Muc4" "Fst"
intersect(row.names(rna_markers[rna_markers$p_val_adj < 0.05, ]),
da_bins_k4[da_bins_k4$p_val < 0.05, ]$gene_name)
#Top 3: "Apoe" "Cxcl14" "Sparc"
roi = c("Cxcl14", "Muc4", "Fst", "Apoe","Sparc")
VlnPlot(seurat, roi, assay = "SCT", group.by = "annot_rna", cols = c(mypal[11], mypal[10]))6 Coverage plots
6.1 Using the epigenomic annotation
DefaultAssay(cre_list[["h3k4me1"]]) <- "bin_50k"
DefaultAssay(cre_list[["h3k27me3"]]) <- "bin_50k"
cre_list[["h3k4me1"]] <- SetIdent(cre_list[["h3k4me1"]], value = cre_list[["h3k4me1"]]@meta.data[["annot_k4"]])
cre_list[["h3k27me3"]] <- SetIdent(cre_list[["h3k27me3"]], value = cre_list[["h3k27me3"]]@meta.data[["annot_k27"]])Sub-setting the max cell for single-cell tracks.
# take one well-covered luminal and basal cell that are non-ambiguously annotated by all annotations
max_cell_k27 <- subset(cre_list[["h3k27me3"]], cells = c("L539C172", "L539C186"))
max_cell_k4 <- subset(cre_list[["h3k4me1"]], cells = c("L539C302", "L539C290")) 6.1.1 Fst
Marker of basal cells by RNA. Different in H3K27me3 in basal-like and luminal-like cells by H3K27me3
#roi="Fst"
#chr13:114,404,115-114,509,113
roi <- GRanges(seqnames = Rle(c("chr13"), c(1)),
ranges = IRanges(114440000, end = 114500000, names = "Fst"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 800) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
p_fst_k27 <- CombineTracks(plotlist = list(cov_plot_k27,gene_plot),
heights = c(15,8)) &
theme(axis.title.y = element_text(size = 7))
p_fst_expression <- VlnPlot(seurat, "Fst", assay = "SCT", group.by = "annot_rna", cols = c(mypal[11], mypal[10]))
p_fst_k27
p_fst_expression
6.1.2 Apoe
Marker of basal cells by RNA. Different in H3K4me1 in basal-like and luminal-like cells by H3K4me1
#roi="Apoe"
#chr7:19,684,744-19,704,653
roi <- GRanges(seqnames = Rle(c("chr7"), c(1)),
ranges = IRanges(19689094, end = 19702731, names = "Apoe"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 800) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
p_apoe_k4 <- CombineTracks(plotlist = list(cov_plot_k4,gene_plot),
heights = c(15,8)) &
theme(axis.title.y = element_text(size = 7))
p_apoe_expression <- VlnPlot(seurat, "Apoe", assay = "SCT", group.by = "annot_rna", cols = c(mypal[11], mypal[10]))
p_apoe_k4
p_apoe_expression
ggsave(file.path(outputDir_plots, "fst_tracks_by_epigenetic_annotation.pdf"),
plot = p_fst_k27,
device = "pdf",
units = "mm",
width = 200,
height = 100)
ggsave(file.path(outputDir_plots, "apoe_tracks_by_epigenetic_annotation.pdf"),
plot = p_apoe_k4,
device = "pdf",
units = "mm",
width = 200,
height = 100)6.2 Using the RNA annotation
cre_list[["h3k4me1"]] <- SetIdent(cre_list[["h3k4me1"]], value = cre_list[["h3k4me1"]]@meta.data[["annot_rna"]])
cre_list[["h3k27me3"]] <- SetIdent(cre_list[["h3k27me3"]], value = cre_list[["h3k27me3"]]@meta.data[["annot_rna"]])
max_cell_k27 <- SetIdent(max_cell_k27, value = max_cell_k27@meta.data$annot_rna)
max_cell_k4 <- SetIdent(max_cell_k4, value = max_cell_k4@meta.data$annot_rna)6.3 Keratines
#Keratine locus
#chr11:100,117,140-100,448,915
roi <- GRanges(seqnames = Rle(c("chr11"), c(1)),
ranges = IRanges(start = 100117140, end = 100339961, names = "Krt"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 4000) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k27_max <- CoveragePlot(
object = max_cell_k27,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 4000) +
ggtitle("H3K27me3 single cell") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 4000) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_k4_max <- CoveragePlot(
object = max_cell_k4,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 4000) +
ggtitle("H3K4me1 single cell") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_rna <- BigwigTrack(
bigwig = list_rna_bw,
region = roi,
smooth = 1000,
type = "coverage",
y_label = "Normalised signal") +
scale_fill_manual(values = rep(mypal[1], 5)) +
theme(legend.position = "none") +
ggtitle("RNA pseudobulk")
p_krt <- CombineTracks(plotlist = list(cov_plot_rna,
cov_plot_k27_max, cov_plot_k27,
cov_plot_k4_max, cov_plot_k4, gene_plot),
heights = c(40,15,15,15,15,8)) &
theme(axis.title.y = element_text(size = 7))
tile_plot_k27 <- TilePlot(
object = cre_list[["h3k27me3"]],
region = roi,
tile.cells = 70,
tile.size = 2000) + scale_fill_gradient(low = "white", high = "#00806c")
tile_plot_k4 <- TilePlot(
object = cre_list[["h3k4me1"]],
region = roi,
tile.cells = 70,
tile.size = 2000) + scale_fill_gradient(low = "white", high = "#1c6779")
tiles_krt <- CombineTracks(plotlist = list(tile_plot_k27, tile_plot_k4),
heights = c(15,15)) & theme(axis.title.y = element_text(size = 7))
p_krt
tiles_krt
ggsave(file.path(outputDir_plots, "krt_tracks_by_rna.pdf"),
plot = p_krt,
device = "pdf",
units = "mm",
width = 200,
height = 240)
ggsave(file.path(outputDir_plots, "krt_tiles_by_rna.pdf"),
plot = tiles_krt,
device = "pdf",
units = "mm",
width = 200,
height = 150)6.3.0.1 Krt19
#roi="Krt19"
#chr11:100,138,496-100,148,799
roi <- GRanges(seqnames = Rle(c("chr11"), c(1)),
ranges = IRanges(100138496, end = 100148799, names = "Krt19"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k27_max <- CoveragePlot(
object = max_cell_k27,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 single cell") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_k4_max <- CoveragePlot(
object = max_cell_k4,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 single cell") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_rna <- BigwigTrack(
bigwig = list_rna_bw,
region = roi,
smooth = 200,
type = "coverage",
y_label = "Normalised signal") +
scale_fill_manual(values = rep(mypal[1], 5)) +
theme(legend.position = "none") +
ggtitle("RNA pseudobulk")
p_krt19 <- CombineTracks(plotlist = list(cov_plot_rna,
cov_plot_k27_max, cov_plot_k27,
cov_plot_k4_max, cov_plot_k4, gene_plot),
heights = c(30,15,15,15,15,8)) &
theme(axis.title.y = element_text(size = 7))
tile_plot_k27 <- TilePlot(
object = cre_list[["h3k27me3"]],
region = roi,
tile.cells = 70,
tile.size = 400) + scale_fill_gradient(low = "white", high = "#00806c")
tile_plot_k4 <- TilePlot(
object = cre_list[["h3k4me1"]],
region = roi,
tile.cells = 70,
tile.size = 400) + scale_fill_gradient(low = "white", high = "#1c6779")
tiles_krt19 <- CombineTracks(plotlist = list(tile_plot_k27, tile_plot_k4),
heights = c(15,15)) & theme(axis.title.y = element_text(size = 7))
p_krt19
tiles_krt19
ggsave(file.path(outputDir_plots, "krt19_tracks_by_rna.pdf"),
plot = p_krt19,
device = "pdf",
units = "mm",
width = 200,
height = 240)
ggsave(file.path(outputDir_plots, "krt19_tiles_by_rna.pdf"),
plot = tiles_krt19,
device = "pdf",
units = "mm",
width = 200,
height = 150)6.3.0.2 Krt14
#roi="Krt14"
#chr11:100,198,634-100,211,215
max_cell_k27 <- subset(cre_list[["h3k27me3"]], cells = c("L539C172", "L539C186"))
max_cell_k4 <- subset(cre_list[["h3k4me1"]], cells = c("L539C302", "L539C290"))
roi <- GRanges(seqnames = Rle(c("chr11"), c(1)),
ranges = IRanges(100198634, end = 100211215, names = "Krt14"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k27_max <- CoveragePlot(
object = max_cell_k27,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 single cell") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_k4_max <- CoveragePlot(
object = max_cell_k4,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 single cell") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_rna <- BigwigTrack(
bigwig = list_rna_bw,
region = roi,
smooth = 200,
type = "coverage",
y_label = "Normalised signal") +
scale_fill_manual(values = rep(mypal[1], 5)) +
theme(legend.position = "none") +
ggtitle("RNA pseudobulk")
p_krt14 <- CombineTracks(plotlist = list(cov_plot_rna,
cov_plot_k27_max, cov_plot_k27,
cov_plot_k4_max, cov_plot_k4, gene_plot),
heights = c(30,15,15,15,15,8)) &
theme(axis.title.y = element_text(size = 7))
tile_plot_k27 <- TilePlot(
object = cre_list[["h3k27me3"]],
region = roi,
tile.cells = 70,
tile.size = 200) + scale_fill_gradient(low = "white", high = "#00806c")
tile_plot_k4 <- TilePlot(
object = cre_list[["h3k4me1"]],
region = roi,
tile.cells = 70,
tile.size = 200) + scale_fill_gradient(low = "white", high = "#1c6779")
tiles_krt14 <- CombineTracks(plotlist = list(tile_plot_k27, tile_plot_k4),
heights = c(15,15)) & theme(axis.title.y = element_text(size = 7))
tiles_krt14
p_krt14
ggsave(file.path(outputDir_plots, "krt14_tracks_by_rna.pdf"),
plot = p_krt14,
device = "pdf",
units = "mm",
width = 200,
height = 240)
ggsave(file.path(outputDir_plots, "krt14_tiles_by_rna.pdf"),
plot = tiles_krt14,
device = "pdf",
units = "mm",
width = 200,
height = 150)6.3.0.3 Krt17
#roi="Krt17"
#chr11:100,252,878-100,265,917
roi <- GRanges(seqnames = Rle(c("chr11"), c(1)),
ranges = IRanges(100252878, end = 100265917, names = "Krt17"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k27_max <- CoveragePlot(
object = max_cell_k27,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K27me3 single cell") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_k4_max <- CoveragePlot(
object = max_cell_k4,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 400) +
ggtitle("H3K4me1 single cell") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_rna <- BigwigTrack(
bigwig = list_rna_bw,
region = roi,
smooth = 200,
type = "coverage",
y_label = "Normalised signal") +
scale_fill_manual(values = rep(mypal[1], 5)) +
theme(legend.position = "none") +
ggtitle("RNA pseudobulk")
p_krt17 <- CombineTracks(plotlist = list(cov_plot_rna,
cov_plot_k27_max, cov_plot_k27,
cov_plot_k4_max, cov_plot_k4, gene_plot),
heights = c(30,15,15,15,15,8)) &
theme(axis.title.y = element_text(size = 7))
tile_plot_k27 <- TilePlot(
object = cre_list[["h3k27me3"]],
region = roi,
tile.cells = 70,
tile.size = 200) + scale_fill_gradient(low = "white", high = "#00806c")
tile_plot_k4 <- TilePlot(
object = cre_list[["h3k4me1"]],
region = roi,
tile.cells = 70,
tile.size = 200) + scale_fill_gradient(low = "white", high = "#1c6779")
tiles_krt17 <- CombineTracks(plotlist = list(tile_plot_k27, tile_plot_k4),
heights = c(15,15)) & theme(axis.title.y = element_text(size = 7))
p_krt17
tiles_krt17
ggsave(file.path(outputDir_plots, "krt17_tracks_by_rna.pdf"),
plot = p_krt17,
device = "pdf",
units = "mm",
width = 200,
height = 240)
ggsave(file.path(outputDir_plots, "krt17_tiles_by_rna.pdf"),
plot = tiles_krt17,
device = "pdf",
units = "mm",
width = 200,
height = 150)6.4 A region with differences in both K27 & K4
roi <- GRanges(seqnames = Rle(c("chr2"), c(1)),
ranges = IRanges(33005895, end = 34356499, names = "roi"))
gene_plot <- AnnotationPlot(
object = cre_list[["h3k4me1"]],
region = roi)
cov_plot_k27 <- CoveragePlot(
object = cre_list[["h3k27me3"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 10000) +
ggtitle("H3K27me3 pseudobulk") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k27_max <- CoveragePlot(
object = max_cell_k27,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 10000) +
ggtitle("H3K27me3 single cell") +
scale_fill_manual(values = c(mypal[3], mypal[3], mypal[3]))
cov_plot_k4 <- CoveragePlot(
object = cre_list[["h3k4me1"]],
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 10000) +
ggtitle("H3K4me1 pseudobulk") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_k4_max <- CoveragePlot(
object = max_cell_k4,
region = roi,
annotation = FALSE,
peaks = FALSE,
window = 10000) +
ggtitle("H3K4me1 single cell") +
scale_fill_manual(values = c(mypal[2], mypal[2], mypal[2]))
cov_plot_rna <- BigwigTrack(
bigwig = list_rna_bw,
region = roi,
smooth = 10000,
type = "coverage",
y_label = "Normalised signal") +
scale_fill_manual(values = rep(mypal[1], 5)) +
theme(legend.position = "none") +
ggtitle("RNA pseudobulk")
p <- CombineTracks(plotlist = list(cov_plot_rna,
cov_plot_k27_max, cov_plot_k27,
cov_plot_k4_max, cov_plot_k4, gene_plot),
heights = c(30,15,15,15,15,8)) &
theme(axis.title.y = element_text(size = 7))
tile_plot_k27 <- TilePlot(
object = cre_list[["h3k27me3"]],
region = roi,
tile.cells = 70,
tile.size = 5000) + scale_fill_gradient(low = "white", high = "#00806c")
tile_plot_k4 <- TilePlot(
object = cre_list[["h3k4me1"]],
region = roi,
tile.cells = 70,
tile.size = 5000) + scale_fill_gradient(low = "white", high = "#1c6779")
tiles <- CombineTracks(plotlist = list(tile_plot_k27, tile_plot_k4),
heights = c(15,15)) & theme(axis.title.y = element_text(size = 7))
tiles
p
ggsave(file.path(outputDir_plots, "region1_tracks_window10000_by_rna.pdf"),
plot = p,
device = "pdf",
units = "mm",
width = 200,
height = 240)
ggsave(file.path(outputDir_plots, "region1_tiles_window10000_by_rna.pdf"),
plot = tiles,
device = "pdf",
units = "mm",
width = 200,
height = 150)6.5 Expression of the genes in these regions
roi_krt = c("Krt15", "Krt19", "Krt14", "Krt16", "Krt17", "Eif1")
p_genes_rna_krt <- VlnPlot(seurat, roi_krt, assay = "SCT", group.by = "annot_rna", cols = c(mypal[11], mypal[10]))
roi_region1 = c("Garnl3", "Ralgps1", "Zbtb34", "Lmx1b", "Mvb12b", "Gm13403", "C230014O12Rik", "Angptl2", "Zbtb43", "Gm13530", "9430024E24Rik", "Gm13404", "Pbx3")
p_genes_rna_region1 <- VlnPlot(seurat, roi_region1, assay = "SCT", group.by = "annot_rna", cols = c(mypal[11], mypal[10]))
p_genes_rna_krt
p_genes_rna_region1
ggsave(file.path(outputDir_plots, "genes_rna_region1.pdf"),
plot = p_genes_rna_region1,
device = "pdf",
units = "mm",
width = 250,
height = 250)
ggsave(file.path(outputDir_plots, "genes_rna_krt.pdf"),
plot = p_genes_rna_krt,
device = "pdf",
units = "mm",
width = 200,
height = 200)